%% load calibration North at initial steady state and solve for the steady state equilibrium in the Northern region
str_region="North_";
str1="params_";
 str_3="calibration_Table2";
 str_save=append(path_calibration,str1,str_region,str_3,".mat");
filename=cell2mat(str_save);
load(filename)
A=A_n_ben;read_pars;read_pars_n;pars_n=pars;

Z_n=z_0_n*reshape(exp(z_vec)'*ones(1,N_b),Ntot,1);
A_n_new=A;
pp_inve_n=pp_inv;
b_common = B0_n/(z_0*A_n_new);
b0_low_n=b_common/zhig_n;
b0_mid_n=b_common/zmid_n;
b0_hig_n=b_common/zlow_n;
b_bar_n=b_bar;
% find invariant distribution (b,z)
b_up=b0_hig_n;b_md=b0_mid_n;b_dn=b0_low_n;
q_hig=q_hig_n;q_low=q_low_n;q_mid=q_mid_n;
bhig=b0_hig_n;bmid=b0_mid_n;blow=b0_low_n;
incid=1;
get_distribution_b;
inde_med = find(cumsum(f_dist)/sum(f_dist)>0.5,1)-1;
med_debt_coupon_n = coupon./ppval(pp_x,b_vec_d(inde_med));

b_vec_mat_l_n = b_vec_mat_l+log(b_bar);
b_vec_mat_n   = exp(b_vec_mat_l_n);
b_mat_n       = ones(length(z_vec),1)*(b_vec_mat_n);
bb_vec_mat_n  = reshape(b_mat_n,N_z*N_b,1);
b_n           = reshape(ones(N_z,1)*(b_vec_mat_n),Ntot,1);
inv_n         = ppval(pp_inve_n,b_n);
b_vec_mat=b_vec_mat_n;bb_vec_mat=bb_vec_mat_n;
policy_expe=0;
get_distribution_bZ;
f_dist_n=(f_dist);
P_default_avg_n   = 1/sum(f_dist)+delta_n;
debt_n            = (bb_vec_mat'*f_dist_n)/sum(f_dist);
z_n               = (Z_n'*f_dist_n)/sum(f_dist);
A_T_n   = A_T;
h_00_n  = h_00;
inde_b_bar_n=inde_b_bar;
inde_bi_n=inde_bi;
inde_bb_n=inde_bb;
inde_nego_n = inde_nego;
A_T_n_1=A_T;
h_00_n_1  = h_00;
sum_n_1=sum(f_dist);

%avg rate paid on bonds
avg_debt_coupon_n = ((coupon./ppval(pp_x,bb_vec_mat))'*f_dist_n)/sum(f_dist);
calib=1;a= (a_targ-(1-phi)*varphi)/phi;
get_firm_policy;

%evaluate condition on equilibrium existence and non-monotonicity of b*x_0(b)
b_max=fmincon(@(b) -b.*ppval(pp_x,b),3)
rat_tresh=ppval(pp_v,b_max)/(ppval(pp_x,b_max)*b_max)
rat_need =(varzeta/(r+delta))/(k-varzeta/(r+delta))

b_vec_n = 1:0.01:4;Ex0_vec=[];
for ii=1:length(b_vec_n)
    b0_mid_n=b_vec_n(ii)/zmid_n;
    b0_hig_n=b_vec_n(ii)/zlow_n;
    x0_hig = ppval(pp_x,b0_hig_n);
    x0_mid = ppval(pp_x,b0_mid_n);
    Ex0_vec(ii)    = min(1,(x0_hig*q_hig+x0_mid*q_mid));
end

figure(300)
plot(b_vec_n*(z_0*A_n_new),b_vec_n.*Ex0_vec)
xline(b_common*(z_0*A_n_new))
